In [1]:
from jedi import jedi
from jedi.utils import plot, seedutil, func_generator, init_tools
from functools import partial
import random
import types
import matplotlib.pylab as plt
from matplotlib import animation
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from __future__ import division
from scipy.integrate import odeint, ode
from numpy import zeros,ones,eye,tanh,dot,outer,sqrt,linspace, \
cos,pi,hstack,zeros_like,abs,repeat
from numpy.random import uniform,normal,choice
import numpy.linalg as la
from ipywidgets import interact, fixed
from sklearn.decomposition import PCA
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
In [2]:
reload(jedi)
reload(seedutil)
Out[2]:
In [3]:
# Setting Seeds
seeds = random.sample(range(100000), 1)
In [4]:
t = linspace(-2,2,1000)
plt.subplot(2,1,1)
plt.plot(t, jedi.step_decode(t));
plt.ylim(-0.5,1.5)
delta = 100
plt.subplot(2,1,2)
plt.plot(t, jedi.sigmoid(delta, t));
plt.ylim(-.5,1.5);
In [5]:
reload(jedi)
reload(seedutil)
Out[5]:
In [6]:
# sine-wave target
target = lambda t0: cos(2 * pi * t0/.5)
In [7]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstart = 0 # learning start time
tstop = 7 # learning stop time
rho = 1.5 # spectral radius of J
N = 300 # size of stochastic pool
lr = 1.0 # learning rate
pE = .8 # percent excitatory
sparsity = (.1,1,1) # sparsity
In [8]:
noise_mat = np.array([np.random.normal(0,0.3,N) for i in range(1002)])
In [9]:
errors = []
xs = []
for seedling in seeds:
J, Wz, _, x0, u, w = init_tools.set_simulation_parameters(seedling, N, 1, pE=pE, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, params):
i = params['index']
tanh_x = params['tanh_x']
z = params['z']
noise = params['noise'][i]
return (-x + dot(J, tanh_x) + Wz*z + noise)/dt
x, t, z, _, wu,_ = jedi.force(target, model, lr, dt, tmax, tstart, tstop, x0, w, noise=noise_mat)
xs.append(x)
error = np.abs(z-target(t))
errors.append(error)
errors = np.array(errors)
In [10]:
T = 300
plt.figure(figsize=(12,3))
ax = plt.subplot(111)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
for i in range(20):
ax.plot(t[:T], x[:T, i], "#696969")
#plt.xlabel('time', fontweight='bold', fontsize=16)
In [11]:
plt.figure(figsize=(12,3))
plot.signal_error(errors, t, tstart, tstop, title= "FORCE (Sin Wave)", burn_in=5)
In [12]:
plt.figure(figsize=(12,4))
plot.target_vs_output_plus_error(t, z, wu, target, offset=1, log=False)
In [13]:
pca = PCA(n_components=2)
pca.fit(x)
pca_x = pca.transform(x).T
In [14]:
reload(plot)
Out[14]:
In [15]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_2dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));2
Out[15]:
In [16]:
np.save("../data/random/pca_x.npy", pca_x)
In [17]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
In [18]:
tmax = t[-1]
tmin = t[0]
t = t[:-1]
def view_rate(time, xs, t, sim):
ti = np.argmax(tv >= time)
plt.title("t=%.2fms" % ti)
plt.plot(xs[sim][ti, :], np.tanh(xs[sim][ti, :]), 'ro',
alpha=.1)
In [19]:
interact(view_rate, time=(tmin, tmax, .1), xs=fixed(xs),
t=fixed(t), sim=fixed(0));
In [20]:
np.save("../data/random/force_sin_xs.npy", np.array(xs))
In [21]:
reload(jedi)
reload(seedutil)
Out[21]:
In [22]:
jedi.dforce?
In [23]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstart = 0
tstop = 5 # learning stop time
rho = 1.5 # spectral radius of J
N = 300 # size of stochastic pool
lr = 1.0 # learning rate
pE = .8 # percent excitatory
sparsity = (.1,1,1) # sparsity
In [28]:
derrors = []
for seed in seeds:
J, Wz, _, x0, u, w = init_tools.set_simulation_parameters(seedling, N, 1, pE=pE, p=sparsity, rho=rho)
def model(t0, x, params):
index = params['index']
tanh_x = params['tanh_x']
z = params['z']
noise = params['noise'][index]
return (-x + dot(J, tanh_x) + Wz*z + noise)/dt
x,t,z,w_learn,wu,_ = jedi.dforce(jedi.step_decode, targets, model, lr, dt, tmax, tstart, tstop, x0, w,
inputs=inputs, noise=noise_mat, pE=pE)
derror = np.abs(z-target(t))
derrors.append(derror)
derrors = np.array(derrors)
In [27]:
plt.figure(figsize=(12,3))
ax = plt.subplot(111)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
for i in range(10):
ax.plot(t[:T], x[:T, i], "#696969")
In [25]:
plt.figure(figsize=(12,5))
plot.signal_error(derrors, t, tstart, tstop, title="DFORCE (Sin Wave)", burn_in=5)
plot.target_vs_output_plus_error(t, z, wu, target, log=False)
In [ ]:
pca = PCA(n_components=2)
pca.fit(x)
In [ ]:
pca_x = pca.transform(x).T
In [ ]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_2dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [ ]:
plt.figure(figsize=(7,3))
interact(view_rate, time=(tmin, tmax, .1), xs=fixed(xs),
t=fixed(t), sim=fixed(0));
In [ ]:
edi.sigmoid(1, np.array(xs))
In [ ]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
plt.title("SFORCE (Sin Wave)");
In [ ]:
In [ ]:
In [ ]:
plt.figure(figsize=(12,3))
plot.cross_signal_error(errors, derrors, t, tstart, tstop,
title="FORCE vs SFORCE (Sin Wave))", burn_in=100)
In [ ]:
reload(func_generator)
In [ ]:
targets = np.load("../data/stability/flipflop/targets_tmax10sec.npy")
inputs = np.load("../data/stability/flipflop/inputs_tmax10sec.npy")
In [ ]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstart = 0
tstop = 5 # learning stop time
rho = 1.02 # spectral radius of J
N = 300 # size of stochastic pool
lr = 1.0 # learning rate
pE = .8 # percent excitatory
sparsity = (.1,1,1) # sparsity
I = 1 # input size
In [ ]:
noise_mat = np.array([np.random.normal(0,.3, N) for i in range(1002)])
In [ ]:
internal_noise = True
errors = []
wus = []
zs = []
for seedling in seeds[:1]:
J, Wz, Wi, x0, u, w = init_tools.set_simulation_parameters(seedling, N, I, pE=pE, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, params):
index = params['index']
tanh_x = params['tanh_x']
z = params['z']
inp = params['inputs'][index]
if internal_noise:
noise = params['noise'][index]
else:
noise = 0
return (-x + dot(J, tanh_x) + dot(Wi, inp) + Wz*z + noise)/dt
if internal_noise:
x,t,z,w_learn,wu,_ = jedi.force(targets, model, lr, dt, tmax, tstart, tstop, x0, w,
inputs=inputs, noise=noise_mat)
else:
inputs += noise_mat[:, 0]
x,t,z,w_learn,wu,_ = jedi.force(targets, model, lr, dt, tmax, tstart, tstop, x0, w,
inputs=inputs)
zs.append(z)
wus.append(wu)
error = np.abs(z-np.array(targets))
errors.append(error)
errors = np.array(errors)
In [ ]:
In [ ]:
int(tmax/dt+2)
In [ ]:
plt.figure(figsize=(12,5))
plot.signal_error(errors, t, tstart, tstop, title= "FORCE (Flip-Flop)", burn_in=5)
In [ ]:
ind = 0
In [ ]:
plt.figure(figsize=(12,5))
if ind < len(seeds):
print("Seed: %d" % ind)
plot.target_vs_output_plus_error(t, zs[ind], wus[ind], targets, offset=0, log=False)
ind+=1
In [ ]:
T = 500
plt.figure(figsize=(12,3))
ax = plt.subplot(111)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
for i in range(10):
ax.plot(t[:T], x[:T, i], "#696969")
#plt.xlabel('time', fontweight='bold', fontsize=16)
In [ ]:
plt.figure(figsize=(12,5))
plt.legend()
plt.subplot(2,1,2)
for i in range(10):
plt.plot(t[:], x[:, i]);
plt.subplot(2,1,1)
plt.plot(t, np.array(z), lw=2, label="output");
In [ ]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
plt.title("FORCE (Flip Flop)");
In [ ]:
pca = PCA(n_components=2)
pca.fit(x)
pca_x = pca.transform(x).T
In [ ]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_2dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [ ]:
np.save("../data/random/pca_x_ff_noise.npy", pca_x)
In [ ]:
"({0}, {1})".format(.1,.2)
In [ ]:
pca_x.shape
In [ ]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstop = 5 # learning stop time
rho = 1.02 # spectral radius of J
N = 300 # size of stochastic pool
lr = 1.0 # learning rate
pE = .8 # percent excitatory
sparsity = (.1,1,1) # sparsity
I = 1 # input size
In [ ]:
internal_noise = True
derrors = []
wus = []
zs = []
for seedling in seeds[:2]:
J, Wz, Wi, x0, u, w = init_tools.set_simulation_parameters(seedling, N, I, pE=pE, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, params):
index = params['index']
tanh_x = params['tanh_x']
z = params['z']
inp = params['inputs'][index]
if internal_noise:
noise = params['noise'][index]
else:
noise = 0
return (-x + dot(J, tanh_x) + dot(Wi, inp) + Wz*z + noise)/dt
if internal_noise:
x,t,z,w_learn,wu,_ = jedi.dforce(jedi.step_decode, targets, model, lr, dt, tmax, tstart, tstop, x0, w,
inputs=inputs, noise=noise_mat, pE=pE)
else:
inputs += noise_mat[:, 0]
x,t,z,w_learn,wu,_ = jedi.dforce(jedi.step_decode, targets, model, lr, dt, tmax, tstart, tstop, x0, w,
inputs=inputs, pE=pE)
zs.append(z)
wus.append(wu)
derror = np.abs(z-np.array(targets))
derrors.append(derror)
derrors = np.array(derrors)
In [ ]:
T = 500
plt.figure(figsize=(12,3))
ax = plt.subplot(111)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
for i in range(10):
ax.plot(t[:T], x[:T, i], '#696969')
#plt.xlabel('time', fontweight='bold', fontsize=16)
In [ ]:
plt.figure(figsize=(12,5))
plot.signal_error(derrors, t, tstart, tstop, title="SFORCE (Flip-Flop)", burn_in=5)
In [ ]:
# Setting seed index
ind = 0
In [ ]:
plt.figure(figsize=(12,5))
if internal_noise:
plot.target_vs_output_plus_error(t, zs[ind], wus[ind], targets, offset=0, log=False)
else:
plot.target_vs_output_plus_error(t, zs[ind], wus[ind], targets+noise_mat[:,0], offset=0, log=False)
In [ ]:
plt.figure(figsize=(12,4))
plt.legend()
plt.subplot(2,1,2)
for i in range(20):
plt.plot(t[:], x[:, i]);
plt.subplot(2,1,1)
plt.plot(t, np.array(z), lw=2, label="output");
In [ ]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
plt.title("DFORCE (Flip Flop)");
In [ ]:
pca = PCA(n_components=3)
pca.fit(x)
pca_x = pca.transform(x).T
In [ ]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_2dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [ ]:
interact(plot.visualize_3dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [ ]:
import seaborn
In [ ]:
In [ ]:
# cross mean signal error
plt.figure(figsize=(12,5))
plot.cross_signal_error(errors, derrors, t, tstart, tstop,
title="FORCE vs SFORCE (Flip-Flop)", burn_in=5)
In [ ]:
# Setting seed index
ind = 0
In [ ]:
errors[0]
In [ ]:
plt.figure(figsize=(12,5))
print("Seed: %d" % ind)
plot.cross_signal_error(errors[ind], derrors[ind], t, tstart, tstop,
title="FORCE vs SFORCE, (Flip-Flop)",
burn_in=5, mean=False)
ind+=1
In [ ]:
# Parameters specified by Abbott 2009.
def lorentz((x, y, z), t0, sigma=10., beta=8./3, rho=28.0):
return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]
In [ ]:
break_in = 500
T = 1501 # period
x0 = np.random.randn(3) # starting vector
t_= np.linspace(0, 60, T)
lorenz = odeint(lorentz, x0, t_)/20
targets = lorenz[break_in:,0]
t_in = t_[break_in:]
In [ ]:
# Visualizing Lorentz attractor
plt.figure(figsize=(12,3))
plt.plot(t_in, targets);
plt.xlim([min(t_in), max(t_in)])
plt.ylim([-1.5,1.5])
plt.xlabel('time', fontsize=14);
plt.ylabel('y', fontsize=14);
In [ ]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstop = 8 # learning stop time
rho = 1.2 # spectral radius of J
N = 1000 # size of stochastic pool
lr = 1.0 # learning rate
pE = None # percent excitatory
sparsity = (.1,1,1) # sparsity
In [ ]:
errors = []
wus = []
zs = []
for seedling in seeds:
J, Wz, _, x0, u, w = init_tools.set_simulation_parameters(seedling, N, 1, pE=pE, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, tanh_x, inp, z):
return (-x + dot(J, tanh_x) + Wz*z)/dt
x, t, z, _, wu,_ = jedi.force(targets, model, lr, dt, tmax, tstart, tstop, x0, w, 0)
zs.append(z)
wus.append(wu)
error = np.abs(z[1:]-np.array(targets))
errors.append(error)
errors = np.array(errors)
In [ ]:
plt.figure(figsize=(12,7))
ind = 0
tstart = 1
plot.target_vs_output_plus_error(t[tstart:], zs[ind][tstart:],
wus[ind][tstart:], targets[tstart:], offset=1, ylim=[-1,2])
In [ ]:
# Visualizing activities of first 20 neurons
T = 500
K = 20
plt.figure(figsize=(12,4))
plt.subplot(211)
plt.title("Neuron Dynamics");
for i in range(K):
plt.plot(t[:T], x[:T, i]);
plt.subplot(212)
for i in range(K):
plt.plot(t[-T:], x[-T:, i]);
plt.xlim(t[-T], t[-1])
In [ ]:
plt.figure(figsize=(12,5))
plot.signal_error(errors, t[1:], tstop, title= "FORCE", burn_in=5)
In [ ]:
ind = 0
In [ ]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
plt.title("FORCE (Flip Flop)");
In [ ]:
pca = PCA(n_components=3)
pca.fit(x)
pca_x = pca.transform(x).T
In [ ]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_3dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [ ]:
In [ ]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstop = 8 # learning stop time
rho = 1.02 # spectral radius of J
N = 1000 # size of stochastic pool
lr = 1.0 # learning rate
pE = None # percent excitatory
sparsity = (.1,1,1) # sparsity
In [ ]:
derrors = []
wus = []
zs = []
for seedling in seeds:
J, Wz, _, x0, u, w = init_tools.set_simulation_parameters(seedling, N, 1, pE=pE, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, tanh_x, inp, z):
return (-x + dot(J, tanh_x) + Wz*z)/dt
x, t, z, _, wu,_ = jedi.dforce(jedi.step_decode, targets, model, lr, dt, tmax, tstart, tstop, x0, w, 0, pE=.8)
zs.append(z)
wus.append(wu)
derror = np.abs(z[1:]-np.array(targets))
derrors.append(derror)
derrors = np.array(derrors)
In [ ]:
plt.figure(figsize=(12,5))
ind = 0
plot.target_vs_output_plus_error(t, zs[ind], wus[ind], targets, offset=1, log=False)
In [ ]:
plt.figure(figsize=(12,4))
plot.signal_error(derrors, t[1:], tstop, title= "SFORCE", burn_in=5)
In [ ]:
plt.figure(figsize=(12,5))
print("Seed: %d" % ind)
plot.signal_error(derrors[ind], t[1:], tstop,
title= "SFORCE", burn_in=5, mean=False)
plot.target_vs_output_plus_error(t, zs[ind], wus[ind], targets, offset=1, log=False)
ind+=1
In [ ]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
plt.title("SFORCE (Flip Flop)");
In [ ]:
pca = PCA(n_components=3)
pca.fit(x)
pca_x = pca.transform(x).T
In [ ]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_3dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [ ]:
In [ ]:
# cross mean signal error
plot.cross_signal_error(errors, derrors, t[1:], tstop,
title="FORCE vs SFORCE (Lorentz)", burn_in=5)
In [ ]: